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ABSTRACT 

Random walks with a general, nonlinear barrier have found recent applications ranging from reion- 
ization topology to refinements in the excursion set theory of halos. Here, we derive the first-crossing 
distribution of random walks with a moving barrier of an arbitrary shape. Such a distribution is shown 
to satisfy an integral equation that can be solved by a simple matrix inversion, without the need for 
Monte Carlo simulations, making this useful for exploring a large parameter space. We discuss exam- 
ples in which common analytic approximations fail, a failure which can be remedied using the method 
described here. 

Subject headings: cosmology: theory - large-scale structure of universe - intergalactic medium - 
galaxies: halos - galaxies; structure - methods: analytical - methods: numerical 

1. INTRODUCTION 

A random walk is a stochastic process consisting of a 
sequence of uncorrelated discrete steps. It has found ap- 
plications in diverse areas ever since Einstein's study of the 
Brownian motion. In cosmology, it has been used mainly 
to model the statistics of halo formation and mergers, in a 
theory sometimes referred to as extended-Press-Schechter, 
or excursion set (see Press & Schechter 1974, Bond et al. 
1991, Lacey & Cole 1993). The idea is quite simple. Take 
a snapshot of the universe at some early time, with some 
given initial density field. Pick a point anywhere in the 
universe. Imagine smoothing the density field around this 
point with a filter, say a top-hat filter, of a successively 
smaller radius. Naturally, (in a hierarchical universe such 
as our own) one expects that the larger the radius of the 
filter, the smaller the mean overdensity would be within 
the filter. This statement is true on average. For any 
given realization, there will be fiuctuations as one varies 
the radius. In other words, a plot of the smoothed density 
versus the size of the filter traces out a random walk, as 
illustrated in Fig. 1. The overall trend is for the smoothed 
density to rise with a smaller filter radius, but with signif- 
icant fiuctuations. 

The excursion set theory postulates a flat threshold or 
barrier (usually motivated by spherical collapse): once the 
random walk crosses the barrier, a halo is declared to have 
formed with a size given by the crossing radius. Here, the 
filter radius R is often denoted by alternative equivalent 
quantities: mass M or variance S. The latter is particu- 
larly useful for Gaussian random initial conditions - note 
that S increases as the radius or mass is decreased. A 
natural question to ask is: given Gaussian random initial 
conditions, what is the probability that a random walk 
first crosses the barrier between S and S + dS, denoted by 
f{S)dSl The quantity f{S) is known as the first-crossing 
distribution. In excursion set theory, the first-crossing dis- 
tribution is directly related to the halo mass function, a 
key quantity of interest. 

Recently, researchers have found it useful to consider 
random walks with a non-flat barrier - following common 
practice, we will refer to this as a moving barrier, in the 
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Fig. 1. — An example of the random walk pattern in the 
excursion set theory. S (or aj^) denotes the variance of 
Sn which is the density fluctuation smoothed on scale R. 
A large S is equivalent to a small R. The dashed curve 
represents a flat barrier, motivated by spherical collapse. 
The dotted curve represents a moving barrier, motivated, 
for instance, by ellipsoidal collapse. 



sense that the barrier 'moves' with S. ^ For instance, 
Sheth et al. (2001, SMTOl hereafter) found that an eUip- 
soidal collapse model suggests a moving barrier for defining 
halos, which produces a first-crossing distribution, or halo 
mass function, that better matches N-body simulations. 
More generally, random walks can be used to model dis- 

^The term moving barrier could be confusing to some readers: 
even in the case of a flat barrier as in the excursion set theory, the 
barrier could move in the sense that it changes with rcdshift. This 
is not the sense in which we use the term. 



1 



2 



tributions of objects other than dark matter halos. For 
instance, recently, Furlanetto et al. (2004, FZH04 here- 
after) used the first-crossing distribution of random walks 
with a moving barrier to model the HII bubble size distri- 
bution in the reionization epoch. 

Existing techniques to analytically compute the first- 
crossing distribution for a fiat, or even linear, barrier can- 
not be easily generalized to a barrier of an arbitrary shape. 
In this paper, we take an alternative approach, in which 
the first-crossing distribution is shown to satisfy a simple 
integral equation for any kind of barrier. This equation 
can be solved in a straightforward manner by inverting a 
triangular matrix. Not having to resort to Monte Carlo 
simulations is useful especially when exploring a large pa- 
rameter space. We also show that the well-known analytic 
flat /linear solutions are reproduced in this approach. This 
is done in §2. In §3, we briefly describe a few astrophysical 
applications, and show under what circumstances (non- 
Monte-Carlo) analytic approximations that are currently 
in use might fail, a failure for which our approach provides 
a remedy. 

2. RANDOM WALKS WITH A MOVING BARRIER 

In this section, we calculate the first-crossing distribu- 
tion of random walks with a moving barrier. The distri- 
bution satisfies an integral equation which can be solved 
analytically only when the barrier is linear. In general, the 
integral equation can be solved by inverting a triangular 
matrix. The resulting solution is shown to agree well with 
Monte Carlo simulations. 

2.1. The First-Crossing Distribution 

In a 1-D random walk model, the probability distri- 
bution of the displacement 5 (which is equivalent to the 
smoothed density in the excursion set theory) is Gaussian. 
We use S to stand for the variance of the displacement. 
The moving barrier is denoted by B{S). The probability 
that the random walk first crosses the barrier at between 
S and S + dS is represented by f{S)dS. Furthermore, we 
define P{S, S)dS as the probability that the random walk 
crosses between 5 and S + dS at S without ever crossing 
the barrier before S (i.e. less than S). Since the random 
walk either crosses the barrier before S or passes {6, S) 
with 6 < B{S), we have: 
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(1) 



Ignoring the barrier, P{S, S) should be equal to Po{S, S) 
which is the normal Gaussian distribution with variance S 
and defined as: 



PoiS,S) 
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(2) 



When there is a barrier, to get P{S, S), one should subtract 
from Pq (5, S) the probability that the random walk crosses 
the barrier at somewhere before S and subsequently passes 
through {S,S). Therefore: 



P{S,S) = Po{5,S)- [' dS'f{S')Po{5-B{S'),S-S') 
Jo 



(3) 



Wc arc now ready to solve for both f{S) and P{6, S) as 
there are two integral equations. Taking the derivative of 
eq. [1] with respect to S, we get: 

Using eq. [3] in eq. [4], we obtain an integral equation for 
the flrst-crossing distribution f{S), the key equation of our 
paper: 

rS 



where 



and 



f{S) = g,{S) + f dS'f{S')g2{S, S') (5) 
Jo 



gi{S) = {^-2^)Po{B{S),S) (6) 



g,{S,S') = [2^-^^^§^]PoiB{S)-BiS'),S-S') 

(7) 

To reach eq. [5], one needs to use the following two re- 
lations: 

[ d5PoiS-B{S'),S-S')]\s->s' = ,: (8) 

J-oo ^ 



and 
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Eq. [8] can be proven using Taylor expansion of B{S) — 
B{S') around S = ^"(assuming B(S) is differentiable) . 
Eq. [9] can be derived by noticing that Po{S,S) satisfies 
the diffusion equation: 



dPp ^ I d^Pp 
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(10) 



Eq. [5] takes a form sometimes known as the Volterra 
integral equation of the second kind. In general, it has a 
unique solution: 



f{S) =g,{S) + y2 dS'g^{S')Gn{S, S') (11 



where 



Gi (5, 5') = ,92(5, S') 



(12) 



Gn+i{S,S')= r g2{S,u)Gn{u,S')dt 

JS' 



2.2. Linear Barrier 

If the barrier B{S) is a linear fmiction of S, g2{S, S') is 
a function only of S — S' . The solution of eq. [5] can be 
written down in closed form using Laplace transformation. 
The Laplace transformation is defined as: 

/■OO 

fit) = / exp{-St)f{S)dS (13) 
Jo 

poo 

hit) = / cxp{- St)g,{S)dS 
Jo 

POO 

Ut) = / exp[-(5 - S')t]g2{S - S')d{S - S') 
Jo 
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Eq. [5] implies that 

Rt) = giit) + Htr92{t) 
Assuming the barrier has the form: 
B{S) =a + bS 



(14) 
(15) 



We have 
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hit) 



(1 



+ 2t 



) exp(-a6 - aVb^ + 2t){16) 



Vb^ + 2t 



Therefore 



f{t) = exp(-a6 - a-\/62 + 2t) 



(17) 



Inverting eq. [17] gives the well-known Inverse Gaussian 
distribution {e.g. Sheth 1998, Sheth & Tormen 2002): 
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(18) 



This confirms that the integral equation [5] does yield the 
correct solution in the special case of a linear barrier. 

2.3. General Cases 

For a general moving barrier, eq. [5] can be solved nu- 
merically on a mesh with equal spacing (see Press et al. 
1992 for more details): 

Si = ix AS, i = 0,1,..., N,AS = ^ (19) 

The integral equation is effectively a set of linear algebraic 
equations: 

AC ■'^^ AC 

fiSi) = gi{S,) + ^ E 92{Si, Sj - -^)[f{Sj) + f{Sj.,)] 



If we treat f{Si) as a vector, we have: 
F = G + MF 



(20) 



(21) 



where Fj = f{Si) and Gj = gi{Si), both of which are N+1 
dimensional vectors. M is an {N + 1) x {N + 1) matrix 
and defined as: 





A, 
A 

Ai,i 




if j > i; 

j if < j < i; 
if J = & z ^ 0; 
if i = j = 0. 



where Aij = ^g2iSi, Sj — We immediately have 

F = (I — M)~-'-G as the solution of the integral equation. 
I is the identity matrix. Since M is a triangular matrix, 
the equation can be solved iteratively in a straightforward 
way: 



fiSo) 
fiSi) 

f{S^)\^>l 



gi{So) = o 

.gi(^i)[l-Ai,i]-i 



[g^{Si)+Y,fiSj){A,j+\j+^] 



(22) 
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Fig. 2. — The first-crossing distribution for a barrier of 
the form: B{S) = 1 + 0.35 + 0.3S^. The dashed curve 
is from a Monte Carlo simulation; the solid curve is the 
solution of eq. [5] . For comparison, we also show the dotted 
curve which is for a linear barrier of the form: B{S) = 
1 + 0.35. The discrepancy between the simulation and 
our exact solution at very small S is due to the small 
number of random walks that first cross at such an S in 
the simulation. 



There is one complication here: since g2{S,S') ap- 
proaches infinity when S approaches S', one needs to de- 
fine (72(5, 5 — ^) carefully when AS is small. According 



to eq. [7], we know: 

g2{S,S')\s->s' 



(23) 



We re-define g2{S,S- ^) in eq. [22] as: 



AC ^ 

52(5, 5 - — ) = — / g2{S, S')dS' (24) 

^ Js-AS 

To compare our calculation with a Monte Carlo simula- 
tion, the barrier is chosen to be B{S) = 1+ 0.35 -|- 0.35^. 
In Fig. 2, the solid curve is the numerical solution of 
eq. [5] using eq. [22]. For comparison, the dotted curve is 
the first-crossing distribution of a linear barrier of the form 
B{S) = 1-1-0.35 - a nonlinear barrier is sometimes approx- 
imated as a linear one in the literature. The dashed curve 
is from a Monte Carlo simulation. Our integral equation 
approach and the Monte Carlo simulation yield consistent 
results. The former is naturally less noisy. 

3. APPLICATIONS IN COSMOLOGY 

In this section, we discuss the cosmological applications 
of random walks with a moving barrier through two ex- 
amples: the halo mass function in the ellipsoidal collapse 
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model and the HII bubble size distribution during reion- 
ization. We will see that an analytic fitting formula for 
the first-crossing distribution proposed by Sheth & Tor- 
men (2002, ST02 hereafter) is a good approximation when 
the barrier is mildly nonlinear. However, the fitting for- 
mula can significantly differ from the exact solution for a 
more general moving barrier. 

3.1. The Ellipsoidal Collapse Model 

As explained in the introduction, in the usual formula- 
tion of the excursion set theory, the halo mass function 
is given by the first-crossing distribution of random walks 
with a flat barrier whose height is determined by the spher- 
ical collapse model: 



n(m)dm = —f(S)dS 
m 



(25) 



where n{m)dm is the number density of dark matter halos 
with mass between m and m+dm,; f{S) is the first-crossing 
distribution; S is the variance of the mass overdensity at 
scale m; p is the mean mass density. This Press-Schechter 
mass function agrees with numerical simulations reason- 
ably well, but a significant discrepancy is found for small 
halos (see SMTOl for details). 

SMTOl points out that such a discrepancy at small 
scales may result from the fact that the spherical col- 
lapse model is over-simplified. Since in general, peaks 
in a Gaussian random field are ellipsoidal (Doroshkevich 
1970, Bardeen et al. 1986), they argue that the forma- 
tion of halos should be better described by an ellipsoidal 
collapse model, in which a moving barrier should be used 
instead of a flat barrier. The barrier adopted by SMTOl 
is: 

B{S) = ^Ssc[l + I3{au)-"] (26) 

where ;/ = S-^^/S: Sgc is the critical overdensity as in 
the spherical collapse model; the constant parameters are: 
/3 « 0.485, a w 0.75 and a « 0.615. Note that Ssc is a 
function of redshift z. It is found in numerical simulations 
that the ellipsoidal collapse model provides a significant 
improvement over the spherical collapse model. 

To find out the first-crossing distribution of a random 
walk model with a moving barrier of the form eq. [26], 
SMTOl suggests a fitting formula which is based on Monte 
Carlo simulations. ST02 proposes a more general fitting 
formula for moving barriers: 



f{S)dS=\T{S)\exp[ 



B{S)^ dS/S 
2S ^^/2^ 



where 



T{S) = Yl 



n=0 



{~S)" d"B{S) 



(27) 



(28) 



and B{S) is the barrier. It is obvious that eq. [27] gives 
the exact first-crossing distribution for a linear barrier. 

We compare eq. [27] with the exact solution of eq. [5] in 
Fig. 3. Instead of using f{S) to represent the first- crossing 
distribution, it turns out to be convenient to change vari- 
able to u and define: 
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F{iy)du = f{S)dS 



(29) 



Fig. 3. — The first-crossing distribution for a barrier of 
the form given by eq. [26], which is motivated by ellipsoidal 
collapse. The solid curve is the exact solution from solving 
eq. [5]; the dashed curve is the fitting formula of ST02, 
eq.[27]. 



The fact that Sgc is a function of redshift z implies that 

f{S) is implicitly a fimction of z as well. The function F, 
on the other hand, has the virtue of being independent of 
z, as can be seen from either eq. [5] or eq. [27]. Fig. 3 shows 
that the fitting formula does provide a good approximation 
to the first-crossing distribution in this case. The relative 
error starts to grow larger than ten percent when v < 1, 
which corresponds to halo mass less than about IO^^Mq 
at 2 = 0. 

3.2. The HII Bubble Size Distribution During 

Reionization 

The evolution of HII bubbles around the epoch of reion- 
ization has been studied using both semi-analytic mod- 
els (e.g. Barkana 2002, Haiman & Holder 2003, Loeb 
et al. 2005) and numerical simulations {e.g. Gnedin 
2000, Sokasian et al. 2003, Ciardi et al. 2003, Sokasian 
et al. 2004). More recently, FZH04 proposes a model for 
such HII bubbles based on the excursion set theory. It 
assumes that the amount of neutral gas being ionized in 
the bubble is directly proportional to the amount of mass 
in the collapsed objects with a virial temperature larger 
than W^K. In other words, the ionized fraction is propor- 
tional to the collapsed fraction. To form an HII bubble, 
one requires the ionized fraction to be equal to one, which 
corresponds to a critical collapsed fraction. With the ex- 
cursion set model, one can relate the collapsed fraction 
with the mass overdensity of the bubble. Thus the critical 
collapsed fraction actually refers to a critical overdensity 
for the formation of HII bubbles. The random walk model 
used in the calculation of the halo mass function can be 
completely carried over to determine the HII bubble size 
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Fig. 4. — The shapes of the barriers assumed in FZH04 
for different redshifts. The upper set is at ^ = 20, and the 
lower set is at z = 16. The solid lines are the true barrier 
shapes, and the dashed lines are the linear approximations. 

distribution. According to FZH04, the barrier can be writ- 
ten as: 

B{S) = b,,-A{a''{mmin)-S\^ (30) 

where S is a^ijn)^ the variance of the mass overdensity of 
bubbles on scale mass m; A is a constant determined by 
the radiation efficiency; nimin is the mass of dark matter 
halos with a virial temperature of 10^ i^, which is about 
10^ M© at redshift 20; Ssc is the usual critical overdensity 
in the spherical collapse model. Note that, as before, B{S) 
is an implicit function of the redshift z, through Ssc and 
nT'min- The HII bubble size distribution is: 

nHii{m)dm = ^f{S)dS (31) 

where n^jj (m) is the number density of bubbles with mass 
m; p is the moan mass density; f{S) is the first-crossing 
distribution for the barrier in cq. [30]. 

Due to the lack of a non-Monte-Carlo method of treating 
random walks with a moving barrier, FZH04 uses a linear 
fit of the true barrier instead: 

B{S) = Bo + BiS (32) 

where Bq = Ssc - A(j{mmin) and Bi = 2^i^~)- ^^S- ^ 
shows the shapes of both the true barrier and its linear fit 
at redshift z = 20, 16 respectively. 

We calculate the first-crossing distribution in three ways: 
one is to use cq. [5], which gives the exact solution; the 
second way is to use eq. [27], the fitting formula given by 
ST02; we also include the first-crossing distribution for the 
linear barrier defined in eq. [32] . The results are presented 
in Fig. 5. One can see all three methods yield similar 



Fig. 5. — The first-crossing distribution for the barriers 
shown in Fig. 4. The lower set is at 2: = 20; the upper 
set is at z = 16. The dashed curves are for the linear 
barriers. Both the solid and the dotted curves arc for the 
true barriers. The solid lines are the exact solutions; the 
dotted lines are from the fitting formula given by eq. [27] . 

results. The fitting formula of ST02 in fact provides a sur- 
prisingly good approximation to the exact solution, with 
a maximum error of less than 10 percent for the relevant 
range of S. The linear barrier approximation is less accu- 
rate, especially at a large S. 

3.3. Cautionary Remarks on the Fitting Formula 

The discussion above suggests that the fitting formula 
by ST02 works well when the barrier is not too complicated 
(i.e. only mildly nonlinear). However, when the barrier 
has a complex shape, the fitting formula may fail. We 
show a hypothetical example in Fig. 7. The barrier is 
chosen to be B{S) = 1 -|- ^{S — 5)^, which is shown in 
Fig. 6. In this example, the first-crossing distribution 
given by the fitting formula differs significantly from the 
exact solution. 

Non-monotonic barriers such as the one above likely 
have useful applications in cosmological problems. Con- 
sidc;r again for instance the problem of the Hll bubble size 
distribution during reionization. In FZH04, the barrier has 
a simple monotonic form in part because recombination is 
ignored. Recently, Furlanetto & Oh (2005) point out that 
recombination disfavors very large Hll bubbles, with an 
upper limit set roughly by the mean free path of the ion- 
izing photons. Since the mean free path depends on the 
average overdensity of a bubble, and since the rms density 
fluctuation is a function of scale (i.e. S), it is conceivable 
that recombination would modify the barrier shape from 
those shown in Fig. 4 by causing an upturn at small S 
(suppressing large bubbles) , perhaps similar to the one in 
our hypothetical example (Fig. 6). Moreover, one should 
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scope of this paper. It suffices to point out that there are 
situations in which commonly used analytic approxima- 
tions, whether it be a simple linear approximation or the 
formula proposed by ST02, are inadequate. Our approach 
via solving eq. [5] by matrix inversion might provide a 
useful alternative, avoiding the need for Monte Carlo sim- 
ulations which can be time consuming when exploring a 
large parameter space. 

4. SUMMARY 

To summarize, we show that the first-crossing distribu- 
tion of random walks with a general moving barrier satis- 
fies eq. [5], the key equation in this paper. This integral 
equation can be solved by inverting a triangular matrix. A 
simple iterative scheme is presented in eq. [22] . We verify 
the technique by comparing with Monte Carlo simulations. 
We also show that the integral equation nicely yields the 
well-known solution for the linear barrier case. Some cos- 
mological applications are briefly discussed, and we cau- 
tion that existing analytic approximations might fail in 
cases with a complex barrier (such as a non-monotonic 
one), where our approach might prove useful. 

Research for this work is supported in part by the DOE 
DE-FG02-92-ER40699 and the NSF AST-0098437. 



Fic;. 6. A hypothetical example of a complex barrier: 
B(5) = 1 + ^(5-5)2. 
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Fig. 7. — The first-crossing distribution of random walks 
with a moving barrier shown in Fig. 6. The solid curve is 
the exact solution; the dotted curve is the fitting formula 
of ST02. 



keep in mind that ionizing sources associated with different 
mass scales (and therefore different S"s) might have quite 
different radiation efficiencies, leading also to the possibil- 
ity of a complex barrier shape for HII bubble formation. 
Going into details of these models would be beyond the 



